﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.IO;

namespace ConsoleApplication1
{
    public class Program
    {
        static void Main(string[] args)
        {
            runProgram(args);
        }

        public static void runProgram(string[] args)
        {
            byte[] wavByteArray1 = File.ReadAllBytes("..\\..\\..\\..\\Sound\\Voice\\1.wav");
            byte[] wavByteArray2 = File.ReadAllBytes("..\\..\\..\\..\\Sound\\Voice\\2.wav");
            byte[] wavByteArray3 = File.ReadAllBytes("..\\..\\..\\..\\Sound\\Voice\\3.wav");
            byte[] wavByteArray4 = File.ReadAllBytes("..\\..\\..\\..\\Sound\\Voice\\4.wav");
            byte[] wavByteArray5 = File.ReadAllBytes("..\\..\\..\\..\\Sound\\Voice\\5.wav");
            byte[] wavByteArray6 = File.ReadAllBytes("..\\..\\..\\..\\Sound\\Voice\\6.wav");
            byte[] wavByteArray7 = File.ReadAllBytes("..\\..\\..\\..\\Sound\\Voice\\7.wav");
            byte[] wavByteArray8 = File.ReadAllBytes("..\\..\\..\\..\\Sound\\Voice\\8.wav");
            byte[] wavByteArray9 = File.ReadAllBytes("..\\..\\..\\..\\Sound\\Voice\\9.wav");
            byte[] wavByteArray10 = File.ReadAllBytes("..\\..\\..\\..\\Sound\\Voice\\10.wav");

            List<byte[]> arrays = new List<byte[]>() { wavByteArray1, wavByteArray2, wavByteArray3, wavByteArray4, wavByteArray5, wavByteArray6, wavByteArray7, wavByteArray8, wavByteArray9, wavByteArray10 };

            List<double[]> convertedArrays = ConvertMultipleByteArraysToDouble(arrays);

            double[] convertedVoiceSampleOne = ConvertByteArrayToDouble(wavByteArray1);

            double[] meanArray = CreateMeanArray(convertedArrays);

            double[] scaledMeanArray = ScaleArray(meanArray);

            double[] scaledVoiceSampleOne = ScaleArray(convertedVoiceSampleOne);

            double[] sumOfSquares = SumOfSquares(scaledMeanArray, scaledVoiceSampleOne);

            String[] stringArray = new String[meanArray.Length];
            String[] stringArrayForSum = new String[sumOfSquares.Length];

            for (int i = 0; i < meanArray.Length; i++)
            {
                stringArrayForSum[i] = sumOfSquares[i].ToString();
                stringArray[i] = meanArray[i].ToString();
            }

            File.WriteAllLines("..\\..\\..\\..\\Sound\\Voice\\meanPCM.txt", stringArray);
            File.WriteAllLines("..\\..\\..\\..\\Sound\\Voice\\SumOfSquares.txt", stringArrayForSum);

            Console.ReadLine();
        }

        public static List<double[]> ConvertMultipleByteArraysToDouble(List<byte[]> arrays) 
        {
            List<double[]> convertedArrays = new List<double[]>();
            for (int i = 0; i < arrays.Count; i++)
            {
                convertedArrays.Add(ConvertByteArrayToDouble(arrays.ElementAt(i)));
            }
            return convertedArrays;
        }

        public static double[] ConvertByteArrayToDouble(byte[] array)
        {
            byte[] headerlessByteArray = CreateHeaderlessByteArray(array);
            Int16[] convertedArray = CreateNewInt16Array(headerlessByteArray);
            ConvertBytesToInt16(headerlessByteArray, convertedArray);
            return CreateNewDoubleArrayAndCopyValues(convertedArray, convertedArray.Length);
        }

        public static List<double[]> ScaleMultipleArrays(List<double[]> arrays)
        {
            List<double[]> scaledArrays = new List<double[]>();
            for (int i = 0; i < arrays.Count; i++)
            {
                scaledArrays.Add(ScaleArray(arrays.ElementAt(i)));
            }
            return scaledArrays;
        }

        public static double[] ScaleArray(double[] toBeScaled)
        {
            double[] absoluteValue = new double[toBeScaled.Length];
            for (int j = 0; j < toBeScaled.Length; j++)
            {
                absoluteValue[j] = Math.Abs(toBeScaled[j]);
            }

            int maxIndex = Helper.MaxIndex(absoluteValue);
            double scale = absoluteValue[maxIndex];

            double[] scaledArray = new double[toBeScaled.Length];
            for (int i = 0; i < toBeScaled.Length; i++)
            {
                scaledArray[i] = toBeScaled[i] / scale;
            }
            return scaledArray;
        }

        public static double[] CreateMeanArray(List<double[]> arrays)
        {
            double[] shortestArray = arrays.ElementAt(0);
            for (int j = 1; j < arrays.Count; j++)
            {
                if (shortestArray.Length > arrays.ElementAt(j).Length)
                {
                    shortestArray = arrays.ElementAt(j);
                }
            }

            double[] meanArray = new double[shortestArray.Length];

            for (int i = 0; i < shortestArray.Length; i++)
            {
                double number = 0.0;
                for (int j = 0; j < arrays.Count; j++)
                {
                    number += arrays.ElementAt(j)[i];
                }
                meanArray[i] = number / arrays.Count;
            }
            return meanArray;
        }

        public static double[] SumOfSquares(double[] array1, double[] array2)
        {
            double[] sumOfSquaresArray = new double[array1.Length];
            for (int i = 0; i < array1.Length; i++)
            {
                sumOfSquaresArray[i] = Math.Pow(array1[i], 2) - Math.Pow(array2[i], 2);
            }
            return sumOfSquaresArray;
        }

        public static double[] RunFFT(byte[] wavByteArray, string filenamePCM, string filenameFFT)
        {
            GetSampleRate(wavByteArray);
            GetBitsPerSample(wavByteArray);

            byte[] headerlessWavByteArray = CreateHeaderlessByteArray(wavByteArray);

            Int16[] converted = CreateNewInt16Array(headerlessWavByteArray);

            ConvertBytesToInt16(headerlessWavByteArray, converted);

            String[] StringRepresentationOfArray = CreateStringRepresentationOfArray(converted);

            File.WriteAllLines("..\\..\\..\\..\\Sound" + filenamePCM, StringRepresentationOfArray);

            Console.WriteLine("DONE");
            Console.ReadLine();

            uint power;
            int size;
            FineNearestPowerOfTwoOfArray(converted, out power, out size);

            double[] real = CreateNewDoubleArrayAndCopyValues(converted, size);

            double[] imaginary = CreateImaginaryArray(size);

            Console.WriteLine("Beginning of FFT");
            Console.ReadLine();

            FFT2 fft = new FFT2();
            fft.init(power);
            var result = fft.run(real, imaginary);
            double sum = 0.0;
            for (int i = 0; i < result.Length; i++)
                sum += result[i];
            sum /= result.Length;
            Console.WriteLine(sum);
            Console.ReadLine();

            String[] StringArrayForFFTResult = new String[result.Length];
            for (int i = 0; i < result.Length; i++)
                StringArrayForFFTResult[i] = result[i].ToString();

            File.WriteAllLines("..\\..\\..\\..\\Sound" + filenameFFT, StringArrayForFFTResult);
            return result;
        }

        public static String[] CreateStringRepresentationOfArray(Int16[] converted)
        {
            String[] StringRepresentationOfData = new String[converted.Length];
            for (int i = 0; i < converted.Length; i++)
                StringRepresentationOfData[i] = converted[i].ToString();
            return StringRepresentationOfData;
        }

        public static double[] CreateImaginaryArray(int size)
        {
            double[] imaginary = new double[size];
            for (int i = 0; i < imaginary.Length; i++)
                imaginary[i] = 0.0;
            return imaginary;
        }

        public static double[] CreateNewDoubleArrayAndCopyValues(Int16[] converted, int size)
        {
            double[] real = new double[size];
            for (int i = 0; i < converted.Length; i++)
                real[i] = converted[i];
            return real;
        }

        public static void FineNearestPowerOfTwoOfArray(Int16[] converted, out uint power, out int size)
        {
            bool done = false;
            power = 1;
            size = 0;
            while (!done)
            {
                size = (int)Math.Pow(2, power);
                if (size > converted.Length)
                    done = true;
                else
                    power += 1;
            }
        }

        public static void PrintResults(byte[] headerlessWavByteArray, Int16[] converted)
        {
            int j = 0;
            for (int i = 0; i < converted.Length; i++)
            {
                //Console.Write("Bytes: ");
                //for (int l = 0; l < 8; l++)
                //{
                //    Console.Write(headerlessWavByteArray[(j + l)] + " ");
                //}
                //j += 2;
                //Console.WriteLine("Int16: " + converted[i]);
                Console.Write(converted[i]);
                Console.ReadLine();
            }
        }

        public static void ConvertBytesToInt16(byte[] headerlessWavByteArray, Int16[] converted)
        {
            for (int i = 0; i < headerlessWavByteArray.Length; i += 2)
                converted[i / 2] = BitConverter.ToInt16(headerlessWavByteArray, i);
        }

        public static Int16[] CreateNewInt16Array(byte[] headerlessWavByteArray)
        {
            int index = headerlessWavByteArray.Length / 2;
            Int16[] converted = new Int16[index];
            return converted;
        }

        public static byte[] ResizeArrayToMultipleOfTwo(byte[] headerlessWavByteArray)
        {
            int mod = headerlessWavByteArray.Length % 2;
            if (mod != 0)
                Array.Resize<byte>(ref headerlessWavByteArray, headerlessWavByteArray.Length + (2 - mod));
            return headerlessWavByteArray;
        }

        public static byte[] CreateHeaderlessByteArray(byte[] file)
        {
            var temp = new byte[file.Length - 44];

            for (int i = 44; i < file.Length; i++)
            {
                temp[i - 44] = file[i];
            }
            return ResizeArrayToMultipleOfTwo(temp);
        }

        public static void GetBitsPerSample(byte[] file)
        {
            Console.WriteLine("Bits per sample: " + file[34] + " " + file[35]);
            Console.ReadLine();
        }

        public static void GetSampleRate(byte[] file)
        {
            Console.Write("Sample rate: ");
            for (int i = 24; i < 28; i++)
                Console.Write(file[i] + " ");
            Console.ReadLine();
        }
    }

    public static class Helper
    {
        public static int MaxIndex<T>(this IEnumerable<T> sequence)
        where T : IComparable<T>
        {
            int maxIndex = -1;
            T maxValue = default(T); // Immediately overwritten anyway

            int index = 0;
            foreach (T value in sequence)
            {
                if (value.CompareTo(maxValue) > 0 || maxIndex == -1)
                {
                    maxIndex = index;
                    maxValue = value;
                }
                index++;
            }
            return maxIndex;
        }

        public static double[] doubleConverter(byte[] bytes)
        {
            double[] values = new double[bytes.Length / 8];
            for (int i = 0; i < values.Length; i++)
                values[i] = BitConverter.ToDouble(bytes, i * 8);
            return values;
        }
    }

    public class FFT2
    {
        // Element for linked list in which we store the
        // input/output data. We use a linked list because
        // for sequential access it's faster than array index.
        class FFTElement
        {
            public double re = 0.0;		// Real component
            public double im = 0.0;		// Imaginary component
            public FFTElement next;		// Next element in linked list
            public uint revTgt;			// Target position post bit-reversal
        }

        private uint m_logN = 0;		// log2 of FFT size
        private uint m_N = 0;			// FFT size
        private FFTElement[] m_X;		// Vector of linked list elements

        /**
         * Initialize class to perform FFT of specified size.
         *
         * @param	logN	Log2 of FFT length. e.g. for 512 pt FFT, logN = 9.
         */
        public void init(uint logN)
        {
            m_logN = logN;
            m_N = (uint)(1 << (int)m_logN);

            // Allocate elements for linked list of complex numbers.
            m_X = new FFTElement[m_N];
            for (uint k = 0; k < m_N; k++)
                m_X[k] = new FFTElement();

            // Set up "next" pointers.
            for (uint k = 0; k < m_N - 1; k++)
                m_X[k].next = m_X[k + 1];

            // Specify target for bit reversal re-ordering.
            for (uint k = 0; k < m_N; k++)
                m_X[k].revTgt = BitReverse(k, logN);
        }

        /**
         * Performs in-place complex FFT.
         *
         * @param	xRe		Real part of input/output
         * @param	xIm		Imaginary part of input/output
         * @param	inverse	If true, do an inverse FFT
         */
        public double[] run(double[] xRe, double[] xIm, bool inverse = false)
        {
            uint numFlies = m_N >> 1;	// Number of butterflies per sub-FFT
            uint span = m_N >> 1;		// Width of the butterfly
            uint spacing = m_N;			// Distance between start of sub-FFTs
            uint wIndexStep = 1; 		// Increment for twiddle table index

            // Copy data into linked complex number objects
            // If it's an IFFT, we divide by N while we're at it
            FFTElement x = m_X[0];
            uint k = 0;
            double scale = inverse ? 1.0 / m_N : 1.0;
            while (x != null)
            {
                x.re = scale * xRe[k];
                x.im = scale * xIm[k];
                x = x.next;
                k++;
            }

            // For each stage of the FFT
            for (uint stage = 0; stage < m_logN; stage++)
            {
                // Compute a multiplier factor for the "twiddle factors".
                // The twiddle factors are complex unit vectors spaced at
                // regular angular intervals. The angle by which the twiddle
                // factor advances depends on the FFT stage. In many FFT
                // implementations the twiddle factors are cached, but because
                // array lookup is relatively slow in C#, it's just
                // as fast to compute them on the fly.
                double wAngleInc = wIndexStep * 2.0 * Math.PI / m_N;
                if (inverse == false)
                    wAngleInc *= -1;
                double wMulRe = Math.Cos(wAngleInc);
                double wMulIm = Math.Sin(wAngleInc);

                for (uint start = 0; start < m_N; start += spacing)
                {
                    FFTElement xTop = m_X[start];
                    FFTElement xBot = m_X[start + span];

                    double wRe = 1.0;
                    double wIm = 0.0;

                    // For each butterfly in this stage
                    for (uint flyCount = 0; flyCount < numFlies; ++flyCount)
                    {
                        // Get the top & bottom values
                        double xTopRe = xTop.re;
                        double xTopIm = xTop.im;
                        double xBotRe = xBot.re;
                        double xBotIm = xBot.im;

                        // Top branch of butterfly has addition
                        xTop.re = xTopRe + xBotRe;
                        xTop.im = xTopIm + xBotIm;

                        // Bottom branch of butterly has subtraction,
                        // followed by multiplication by twiddle factor
                        xBotRe = xTopRe - xBotRe;
                        xBotIm = xTopIm - xBotIm;
                        xBot.re = xBotRe * wRe - xBotIm * wIm;
                        xBot.im = xBotRe * wIm + xBotIm * wRe;

                        // Advance butterfly to next top & bottom positions
                        xTop = xTop.next;
                        xBot = xBot.next;

                        // Update the twiddle factor, via complex multiply
                        // by unit vector with the appropriate angle
                        // (wRe + j wIm) = (wRe + j wIm) x (wMulRe + j wMulIm)
                        double tRe = wRe;
                        wRe = wRe * wMulRe - wIm * wMulIm;
                        wIm = tRe * wMulIm + wIm * wMulRe;
                    }
                }

                numFlies >>= 1; 	// Divide by 2 by right shift
                span >>= 1;
                spacing >>= 1;
                wIndexStep <<= 1;  	// Multiply by 2 by left shift
            }

            // The algorithm leaves the result in a scrambled order.
            // Unscramble while copying values from the complex
            // linked list elements back to the input/output vectors.
            x = m_X[0];
            while (x != null)
            {
                uint target = x.revTgt;
                xRe[target] = x.re;
                xIm[target] = x.im;
                x = x.next;
            }

            return xRe;
        }

        /**
         * Do bit reversal of specified number of places of an int
         * For example, 1101 bit-reversed is 1011
         *
         * @param	x		Number to be bit-reverse.
         * @param	numBits	Number of bits in the number.
         */
        private uint BitReverse(uint x, uint numBits)
        {
            uint y = 0;
            for (uint i = 0; i < numBits; i++)
            {
                y <<= 1;
                y |= x & 0x0001;
                x >>= 1;
            }
            return y;
        }
    }
}
